General Disclaimer 


One or more of the Following Statements may affect this Document 


• This document has been reproduced from the best copy furnished by the 
organizational source. It is being released in the interest of making available as 
much information as possible. 


• This document may contain data, which exceeds the sheet parameters. It was 
furnished in this condition by the organizational source and is the best copy 
available. 


• This document may contain tone-on-tone or color graphs, charts and/or pictures, 
which have been reproduced in black and white. 


• This document is paginated as submitted by the original source. 


• Portions of this document are not fully legible due to the historical nature of some 
of the material. However, it is the best reproduction available from the original 
submission. 


Produced by the NASA Center for Aerospace Information (CASI) 



551-70-15 

PREPRINT 


NASA I«X- ^3^70 


VARIATION OF PARAMETERS 
; AND THE 

LONG-TERM BEHAVIOR OF 

planetary ORBITERS 

PARTI , , 

THEORY 

■' " ' ' ■' 

BERWARD KAUFMAN 


f 'i 

, . . (. 


GSFC 


FEBRUARY 1970t 








f ‘ f *<f} 'i 


%W:W 


GODDARD SPACE FUGHfGEiFEK 

GREENBELTv MARYLANBf^ 


fr 

•Ojif,', 



c ' -'if.'/' \ ' L-iif, -fU' - 

Ji- ■ , N, 'V iY ,^%w“ 

i'. ^ 


o (ACCESSION NUMBER) ' * 

S _ 

S (pages) 

g' (NASA CR OR TMX OR AD NUMBER) 


(^RU) 


I I ) I 

I y ' " 1 'll 

I i 1 .' J r 

?■ yi" " 



(CATEGORY) 


1 


X-551-70-15 

PREPRINT 


VARIATION OF PARAMETERS AND THE LONG-TERM 
BEHAVIOR OF PLANETARY ORBITERS 


Bernard Kaufman 

Mission and Systems Analysis Branch 
Mission and Trajectory Analysis Division 


February 1970 


GODDARD SPACE FLIGHT CENTER 
Greenbelt, Maryland 


aiHK or 


fo • 






PRECEDING PAGE BLANK NOT 


VARIATIDN OF PARAMETERS AND THE LONG-TERM 
BEHAVIOR OF PLANETARY ORBITERS 

B. Kaufman 

ABSTRACT 

In this volume, the first of two, is contained the theory 
develc^ed to determine the behavior of planetary orbiters 
under the perturbative influences of the oblat^ess of the 
central body, the presence of the Sun and atmospheric drag. 
Tlte significant ccmtribution of this study is the treatment of 
the third body which includes the medium-period as well as 
the long-period terms of the disturbing function in a system 
of singly- averaged equations. The variation of the elemeits 
is then expressed along with the effects of drag and oblate- 
ness in a of equations to yield a time history of the 
orbit. These equations are easily adapted to an electronic 
computer. In test runs this program has been shown to be 
extremely fast and accurate to within short-period terms 
wh(^ compared with a complete numerical integration of the 
equations of motion. 
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VARIATION OF PARAMETERS AND THE LONG-TERM 
BEHAVIOR OF PLANETARY ORBITERS 

By 

B. Kaufman 

Goddard Space Flight Center 


INTRODUCTION 

This investigation is done in two parts. The theory of the orbital evolution appears in this 
document along wiUi several sample cases. The detailed study of a parametric nature appears in 
Volume n. 

The desi/gn of any mission to place a satellite in orbit about a planet must include a detailed 
investigation of the possible types of orbits that would be useful. These orbits will to a large ex- 
tent be determined by the type of scientific investigations to be carried out in the vicinity of the 
planet. Of particular importance will be the distance of closest approach Opoxlapsis) and the fur- 
thest distance reached (apoapsis). If periapsis is very low^ then it is likely that the orbit will be 
affected by the atmosphere as well as the oblateness of the planet. Wh^ the ap^ipsis becomes 
large, then perturbations caused by other bodies can become significant. As one can easily see it 
is possible for the orbital motion to become higdily complex when one ccmsiders the effects of oblate- 
ness, drag and other bodies and the possibilities of coupling betwemi them. 

It is precisely because of this complexity that a time history of the size, shape and orientation 
of the orbit should be obtained in order to g^dn the maximum sci^tific data from the sateUite, l^ch 
a history should include not only lifetime predictions but also a reasonably accurate history of all 
of the orbital parameters for a variety of initial conditions. Because such a variety of conditions 
will be used, it also becomes essential that the model chosen to produce the time history be not 
only accurate but very fast Any good n-body precision integraticm program of the Cowell or Encke 
type is capable of meeting the first of these criteria but certainly not the second. R was to meet 
both of the above criteria that the present investig^on was undertaken. 

THE SOLAR POTENTIAL: 

The first problem to be considered here is that of third-body perturbations, bi the context used 
here this means the influxes of the Sun on a satellite in orbit about a planet. Although this study 
doas not incluc^ the Earth, the theory presented here can be used by inclucbng both the Sun and the 
Moon as disturbing bodies. To determine the exact effects of third-body perturbations, it is necessary 
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to solve the three-body problem. This would normally require numerical integraticm of the equa- 
itions of motion which, as pointed out earlier, becomes very time consuming when api>lied over 
periods of years for each orbit. Previous experience with such integrations has shown tlmt in 
many cases the amplitudes of the short-term variations in the orbital elements are very small 
wh(m compared with the values of the elements themselves. This means that one may "average" 
the disturbing function over the period of the satellite thereby eliminating all short-perioo terms. 
ThiCi averaging may be accon plished by severiU means such as app^r in References 1 and 2. It 
is also possible to average this result once again over one revolution of the central body about the 
disturbing body (the Sun ir our case). This process thus eliminates all "medium"-period terms 
leaving only the long-period perturbations for consideration. This doubly-averaged system was 
used in References 2, 3, and 4. As pointed out in Reference 4, however, this double-averaging tech- 
nique has some severe limitations, especially for Mars where the interactions ot obUlteness and 
tliird-body effects become very complex. Therefore, the model used in this report is one of a set 
of singly-averaged equations which have been developed in a form which makes diem readily adap- 
able to efficient programming. 

The development of the solar disturMng function can be found in Reference 5, which was the 
preliminary investigation that motivated the present more detailed model. Because of the im- 
portance of this function in deriving the variaticmal equations in a most useful form, tiie develop- 
ment will be r^eated here. 

With reference to Figure 1 we may write the solar cUsturfaing function as follows: 



SATfLLiTE ' -*► where 



we may rewrite Equatimi 1 as 
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For the cases that we will be considering^we note that r/r' « i and therefore we may expand the 
square root term in Equation 2 and neglect a^,l terms of order (r/r*)^ and higher. This yields 




1 / 2r \ 

1 3 , 

/4r* 2 \ r 

R ' 

r' 

V 

• "T^— '°«Sj 

1 +-g- 1 



or 


. I 3 ^ „ 

~ "T7^ +T772 s 


r 


(4 l)] 


t - 




(3) 


The motion of any satellite under the sole influence of a third body can be completely described 


by 


V = 


{*) 


where r is the radius vector from the central body to the satellite, R is the disturbing function and 
V is the gradient operator. This gradi^t is with respect to the state of the satellite but the first 
term of Equation C contains only elements of the Sun. It is obvious then that 


-* M® 

7“7 « 0 

r 


and we can drop this term and write 


R ' 
*8 


.ill!! r i1 

2r'3 rr' / M 


or alternatively 


= ITT [3 “•’« -»]■ 


(5) 


Letting n ^ a where n' is the mean motion ot the 8un and <»' is the semimajor axis, we may 
write Equation 5 as 


R 


0 




Now 


cos S 



*0 • 


(7) 


and 


r® = P cos 0 + Q sind^ 


( 8 ) 


where e' is the true anomaly and 



cos Qcos - sinH sirut^cos i 
sin Hcos <*) + cos Q sin cu cos i 
sin i sin CO 


( 9 ) 


Q 


“ cos n sin CO - sin n cos cocos i 
“ s in n sin CO + cos H cos cocos i • 
sin i cos co 


( 10 ) 


p is a unit vector in the satellite's orbit plane from the center of the planet to pericenter, Q is a 
unit vector in the orbit plane perpendicular to P and in the direction of satellite motion, and i is the 
orbital inclination. Tor completeness and because we shall need it later, we define the third unit 
vector R as the cross product of P and Q 


R 


P X Q 


sinH sin i 
cos n sin i 
cos i 


( 11 ) 


It is important to observe that P, Q and R lure independent of 0, Also that 


0 s 


cos n' cos 
sin n* cds 
sin<3S»' 


( 12 ) 
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wh«re n' is ths light asc&nsioti of tho Sun and is the latitude of the Sun measured from the 
equator of tlie planet. Figure 2 shows the geometry inwlved in tlie above equations and the trans** 
formations to the equators of Venus and Mars appear in Appendix A. 
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Figure 2->Pianet-centered geometry. 
Equation 7 may now be written as 

cos S = (p • 7' ^ cos $ (q • 7' ®) sin $ 


and denoting 


we have 


P * 7'® = a 

Q • - fi 




> 




( 18 ) 


cosS = acos^ + 0sin6 (14) 

where once again we point out that a and are indepmident (d d. Substituting this into Equation 6 
we have 


Re' 


"5 [3 cos* & + 2c0 sin $ cos 6+0^ sin* - ij . 
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Observing that 


l + C08 2 $ 
5 


“ co*^ $ 


t “ co» 26 
J— 


«in^ & 


and 


sin 2^ - 2 sin ^ cos ^ 


we may write the terms in the brackets as 


[a (.» i 


+ COS 26 » . _ ^ -I 1 - cos 26 \ , 

+ o^sin25+/S^ ^ — )-l 


” |( T T T 


and after substituting this, we then have 


r; 


^“5 (a) [(f ii “ l) + T 2^ + sin 2^ I . 


(15) 


As explained earlier, Equation 15 will now be ayeraged over tiie period of the satellite, 
thereby eliminating any dep^idmice of the disturbing function <m the mean anomaly M. This not 
only simplifies the equaticms of motimi but also eliminates all short-term variations over the orbital 
period. We consider therefore only those terms in Equation 15 dependent on e and find that we 
need to integrate three terms ovez^ the period. The three integrals are as follows; 




1 

dft 


Jo 



COS 29 cM 


I 

5^ 



tl sin2^dll 


(18) 


(17) 


(18) 


"i 1 i 


In evaluating these integrals, we will make use of the following relationships 

7 - c 1 - c cos E) 

M “ E ~ e sinE 
clM - dE(l -e cosE) 


where E is the eccentric anomaly. 

Using these equations, we may write Equation 16 as 


(a) ^ ^ I (l“ecosE)^dE 




2n 


^ I (1 “3e cosE -e^ cos^ E + 3e^ cos^ E) dE 


1 2n r . . T 

. . /l+cos2E\ 

1 +3e y 5 j dE 


Evaluating ttiis mtegral, we obtain 


(=) 


dM 


1 +T 


a result valid for all c . 


The integiation of Equation 17 is much more complex 


2n ! (r) cos 2^ dll = ^ (-^) [cos* 0 - sin^ 0] dM 

Jo Jo 




(M 


(19) 

( 20 ) 
(21) 


( 22 ) 
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1 /V 1 

ST J (T/ cos2/i'dM = J [(eos E- e)^ - (l - c^) sin^ eJ clM 


[‘ 


^ j cos^ E- 2 ecos E+ “ (l ~ e*) sin^ E 


dM 


2 ^ \ * 

if 1 + COS2E ^ ^ U-e2) 

^ ^ 2ecosE + e^ j— ^ (l-cos2E) 

Je ^ 


dM 


^ J [(2 ^ - 2 e cosE + (|- + 


COB 2 E 


dM 


1 r'" fse^ 

2;; ■y 

Jo L 


s hr “2 ocose + 


(1 - Y^ cos 2Ej 


dM 




l 3 e^ 

2tt I |Y“"2ecosE + 
'0 


(}" ^2 ) 


(1 - e cos E) dE 


Drofpping all periodic terms which will go to zero, we have 



cos 26 dM 



+ 2e^ cos^ E dE 




which yields 





cos 26 dM 


5 e^ 

2 


(as) 


again a result valid for all e. 


8 



For the last integral, Equation 18, we write the integrand in the following manner: 


[2 sine? cos a'] 




cos td 


2 (l sin E cos E - 2 (l - e^) ‘ ® n sin E 


We then have 


(l~e^) t ^ f' 

Jo 


sin E cos EdM - 




1/2 

I si 

Jo 


in EdM 


(l 


2 ) 12 

I (sin Ecos E “ e sin E cos^ E) dE 

Jo 


(l -e^) 


1/2 C‘ 
Jo 


(si.i E“ e sin Ecos E) dE 


and convenieitly all terms go to zero, leaving 




in 2^ dM * 0 . 


We now define the solar disturbing function as 


1 f’" 

*. = s; 

Jq 


and substituting into this the results of evaluating the integrals, we finally obtain 


^ (~.) [(f'"" +T/32-1) (1 +-^)+( 4 a 2 


a result that is considerably simpler to use than was Equation 15. 
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VARIATION OF PARAMETERS 


We shall here merely write down the equations for the variations in the Keplerian elements. 

The derivation of these equations may be found in any textbook on celestial mechanics (Reference 6) 


^ ^ 3JR 

dt ~ na (5M 


(26) 


dt ena^ ena^ d<x> 


m 


dt 


1 (5R 

na"^|l “ )*^^ sin i 


(28) 


dl CSC i . 3R 

■<JT " ""”2(1 IP " ^ dc^ 


(29) 


duj 

dt 


(1 - e^) (5R cos i 

ena^ ^ na^ (l - )^'^* sin i 


(80) 


dM .. _2_ ^ (l -e^) aR 

^ na da " ^^2 de 


(3.1) 


where in the present case R = R<g. 

Equation 26 yields immediately a very important result Stoce we have averaged out all de> 
pendence of R^^ on we have 


da 

It " 


0 


02) 


which means that there are no long-term or secular variatians in the semimajor axis due to the_ 
presence of a third body. 


From Equation 27 we can derive the variation in eccentricity as follows: we first note that 
as above 



and therefore 


de 

dt 



Inf!® 

2 dct> 
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jEleferring to Equation 25 we have 


^ dR^ <?R^ 

MWIMM ^ . IJ.lTJjr.J . il l»|MIW X ' -t w *■' 

aoi Oa do’ dp acu 


# 


therefore 


da 


ap 

“tt' "■ • 

doj 


0 


and from Equation 9 it is easily seen that 


dp 


Q. 


Thus 


Therefore 


da 




yS = Q 


•» 0 


dyS 


ag 


•jf » 0 


= - p 


VO = 


- a 



. dR 


^3) 


Now 



y Gi4) 
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and subiitituting these into Equation 33 and then into the equation for c , we have 


or 


de (l. 3a^n'^ 




nea' 


2 ■** 4e ^ ) - a/3 ^1 ** e j 


(l -e2)*/2 3n*2^aiy 

ne 2 


(5a^2) 


de 

dt 


15 


(l -e^) 


1/2 



(35) 


For Equation 28 for tiie node, we have 


dn 


1 

na^ (l - e^) sin i 


da di 


aRg dfi 
1J di 


(36) 


where all we need now is da/d i and d^'d i. Now 


and 


da 

dJ 



— ♦ 
ap 

IT 


'■ sinfl sino) sin i 
“ cos ri sinoj sin i 
cos i sin co 


an4 using Equation 11 we have 


da 


R 


# 0 • 

r ” sinoj . 


(37) 


In a similar manner we find 


_ 

IT - 


R • 


r ” cos 0} 


(S8) 
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Combining Equations 34, 37 and 38 and substituting into Equation 36, we have 


en 

Tt 


na 


^ {l - e^) ^ ^ Bin i 




r ' sin 03 


-5“ 30 (l - e2) (g • 0) 


cos u’> 


= r [=>(1*46^) ecu.] 


and letting 


y = R • r'® 


(39) 


we have finally 


dt 


,» 2 


n (l “ e^) sin 


. (| |^a'y(l+4e^) sincy +/?y (l - e^) cos cJj 


(40) 


Next for di/dt, we note that the only term not derived is dR^jdQ. 


dR^ dR^ 5Rg, 

~w = T^m^ujn 


(41) 


where 


da dP _ 

in ~ an ‘ 


* 0 


and 


dP 

in - 


sin n cos CO - cos n sin £jcos i 
cos n cos O' - sin n sin cocos i 
0 


P. 
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Therefore 



( 42 ) 


Similarly 


and 


Then 


a/3 B 
■5TT = a 


•( 0 


5Q - 

- 


sinH sin if - cos Hcos wcos i 
cos n sinw - sin ncos wcos i 
0 




(43) 


Substituting Equations 34^ 42 and 43 into Equation 41 and then alcmg with Equation 33 , the equation 
for inclination becomes 
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3 n*^CiC i /a*Y 
2n(l IrW 


Pj * r'® 


+ ^(l-e^)| • r* ® ~ajS(l + 4e^) cp« i + cy9(l “ c^) co« i 


211(1-02)1/2^7 »ini <^(l+4e2) 


p . 

* X * 


+ /3{l-e2) • r'O - Sa^e^co% i . (44) 


For (ko/dt , we need only to find (9K,/<?e. Since a and 0 are indepaident of e, then from Equation 25 


for we have 


a^n'^/a'V 


- 


3e (ya2 +-|-/32 - + 5 e ^ ^2 - |”/32)j 


[l2ea2-3e;S2.3,] 


a 2 n' 2 /aA'* 


[ 3 ® (4a2 - ^2 _ ijj 


and since 9k.jBi Ims already been evaluated, we have 


duj _ (l-e2)i^^ a2n*2/a*\^ 


V 


i 2 n' 2 /a '\3 


X (r • r'®) sinw + 3 — f~7j 3)S(l -e2) (r • r'®) cos co 
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(4a^ - 1) 


d«u> 

cIF 





n 



3 


3 n'^ 

" Tiry 


(l - e^) 


C08 1 


1/2 


sin 1 



ja(l +4e^) sinw + /? (l ~ e^) coso-j 


rinaUy 


dtt' 3 n'^ 

at ^ nr 



( l - e 2)’»2 


(4a^ “0^ - 1) 


cos t 

^ (l sin i 




Fo^* the last equation dM/dt, we can I '^adily calculate the only needed term as 



Substituting Equations 45 and 47 into Equation 31, we have 
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For convonionce the six equatiotis are writhm below: 









where the subscri|)t meaas third-body variations only. 

This completes die develqpment of the variational equatlcms for third-body perturbations only. 


DRAG 


For the variations due to the presence of an atmosphere, it will be oonv^ient to transform 
Equations 26 through 81 to a slightly different form in which the perturbing force may be written 
in the following componeniform 

P * su, + TU^ + wu^ (50) 

‘ i 

where s is the component along th' radius vector u^ ; x is the component along the normal to 
r, in the orbital plane and such it nakes an angle less dian 90° with the velo(dty vector; and w is 
the compon^jt along the pmrpmidi^ula to the cx^blt plane where x . It can be shown from 


Equations 9, 10 and 11 and thn above definitions that and )mve the foilowing representations 


^co8 Qcor u- «in0 sinu cos i 
U. - [ sin Hcos u-<- cosHsinu cos i 
sin u sin i 


(51) 


U, 


cos fi sinu - sin 0COS u cos i\ 
- sin 0 sinu + cos Hcos u cos i 


cos u sin i 


(53) 


U, 


sinH sin i 
- cos H sin i 
cos i 


(53) 


where u - w+ From Equation 4 we see that the partial derivatives aR/ax, aR/ay, aR/az are the 
components of the acceleration due to the disturbing function R. That is 


aR T* aR aR - 
VR 


m the development of the variation, of parameters, it can be shown that if Cj represents any one of 
the six elem^ts then 


aR ar 

■?C7 ^ 


(i = if • • * . 6) 


(54) 


where 


It is obvious tlmt 


and therefore we need only evaluate 


r = xi + yj + sk . 


F = VR 


aR 

iq- 


3 

^ •'ser 


( 55 ) 
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in terms oi s , T and w and then use tliese partials of R in Equations 26 tlirough 3 1. We will derive 
here only »vr Oa and the derivations of the other partials may be found in Reference 7. From 
Equation 8 we may write 


r cos ti* P + r sin 6Q 




( 56 ) 


Thus 


Then Equation 55 yields 


dr 

“5a 



dR 

~5a 


^ di 



(67) 


Now we can write Equation 56 above as 


where 


Then 


also 


Then 


7 = a(cos E~ e)P ^ a (l ~e*) sin EQ 

^ - r cos <9 = a(cosE“e) 'I 

17 = rsin0 = a ( 1 sin E . 


= COB E - e = - 




= JL 

a 


« “ Pcos0 + Qsind 


Va - - Psin0 + Qcosd 




(58) 


(59) 


(60) 


(61) 


F ‘ P = S(U^ - P) + t(U^ - p) + W(U^ -P) 
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I 


or 


F • P - S cos $ - T sin0 . 


(62) 


Similarly 


F ' Q = S(G, • Q) + T(U^ ' + W(0 ^ ’ q) 


or 


■•-4 

F • Q 


S sin 6 + t cos 6 


and from EquatioraS &9 we have 




F • P - 7-(S^~T7 j) 


and 

F ■ Q = "~(Si 7+T^) . 

Substituting Equations 60 and 64 into ^nation 57, we have 


5R 1 , ^ i 1 V 

Ti- = -T(S^-T^)T + T(S^+T^) a 


“ Tirj) + (St 7^ + T^7?) 


ra 


or 


ai 

3 b 


In a similar manner we can dmive the other partials as; 


as 

3 b 


cos a + Ta sin i 


I + 


i(l - a®) 


[ 



(67) 


dU ^ Sea Ta^ (l - e^) * ^ 

= Tr cos 1 - Wr sin i cos (u +£^) 




Tr 


dR 

•37* s Wrsin(ui + tO. 

Substituting Equations 66 through 70 into Equations 26 through 31, we have 

a (1 - e ' 

— T" 


da 


de 

It ^ 


n(l-e“)l 2 (* 


J*) \ 

Se sint^ + — — p — - Tj 


2 r ^ 1 

j Is sin + T(cos E + cos t/)l 


dn _ r sin (u^+ 6/) 

na2(l-e2)^'2 gin i ^ 


di r cos (to+ $) 

na2(l-s2)l.-! * 


( 68 ) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 




do) 

dt 


cos ^(1 - e2)^'^'2 sin^(l-e2)^^ 

a ne ^ a ne 




T 


r sin (cj + 0) cot i 

a 2 n (l - ^ ^^2 


<M 


n + 


(l “ e^) cos 6 

g-j_g 



S 


(l - e^) sin 8 
i"ne 



where E is the eccentric anomaly. 


(75) 

(76) 


The above form of the variational equations is known as Gauss's form, fii the previous equa- 
tions the variations of the elements were expressed in terms of the partial derivatives of the dis- 
turbing fimction with respect to the elements, fii the Gaussian form we have the variations expressed 
in terms of the disturbing acceleration components. We can further resolve the acceleration in a 
different manner which is useful for the inclusimi of atmospheric drag. We introduce components 
of the acceleration defined as follows 


T' component tang^tial to the orbit 


N compon^t normal to T* and positive toward the central body and 
W same as before. 
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As shoTvn in Reference 8, we then have 


1 + e cos tj ^ e sin 6> 

fl + + 2e cos ^ ^ ^ (l + e“ ^ 2e cos ^ 


e sin & 

(l + + 2e cos 6 


1 + c cos 6 

(l + + 2e cos $Y ^ 


In the present study we consider no lift forces to be present and the drag force acting as a 
negative tangential component (opposite to the velocity vector). We then have 


1 + e cos 6 


^1 + + 2e cos 8 ) 


I '1*1 

1/2 * 


e sin^ 


(l + + 2e cos 


W = N • 0 . 


The drag force is the force ^er unit mass m of magnitude F and 


CpApV^ 


= T' 


where v is the speed, 

Cp is the aerodynamic dra^g coefficient, A is the cross-secUoiial area of the satellite and p is 
the density. Utilizing Equations 73, 74 and 78, we see immediately: 

di) di 

•ar = “ar * 0 . (s 


Fur the other equations we obtain after substituting Equations 78 


da 


‘n'(l-e2y/2 (l^-e2+2ecos0)*/a 


de _ _ ^ p (l -e^) (cos g + e) 

na (l + 4 2e cos 6'j ^ 
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(83) 


iiM ' 

Tir 


^ (l - ^ sin 

n a<‘ • 0 ^ + 2e cos d 


m 


ilM 

OF 


n + 


Cp A y2 sint^ r j + i- e cos t> 1 

nao(i tecost^) ^ ^2 +2e cos(^)i 


m 


(84) 


The presence of sin t' in Equations 83 and 84 implies that du/dt and dM/dt are periodic with small 
amplitudes because of the coefficient Cp A^y m and therefore will be neglected in this treatment. 
However a and e change secularly and must be included. We now make use of the fact that 


2 d^ 


r X V 




and obtain 


dt _ (l - ^ 

^ n(l+ecos(9)^ 


(85) 


also using 


r 


a (l - e^) 

1 + e cos B 


and substituting into 


V2 


6^ 


we have 


y2 




5 " (l + + 2e cos d\ . 

1 -e‘‘ ' 


( 86 ) 


Making use of Equations 85 and 86 in Equatians 81 and 82 we have 


da 

'W 


/x(l + e cos 5)* , 


(87) 


de 

■3^ 


^ (1 - e^) p V(e + cos 0) 
(1 + e cos 0)^ 


m 


( 88 ) 
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where 


M ‘ 0^ . 

We now average the above equations over the orbit and letting a and e represent the perturbations 
in a and e over the orbit, we have 


ni/i J ^ j 4 . ^ dt^ 


e cos H)‘ 


(89) 


CpA(l~e2)^2 r' 




m 


r 

a ^ 


cos 
e cos 


dt' . 


(90) 


The integration is done using Gaussian quadratures, which is very fast and reasonably accurate. 
The use of quadratures was suggested by Uphoff in Reference 4. The density p is taken from 
several models and is calculated as a function of altitude. These models can be found in Refer- 
ences 9 and 10. 

It is to be stressed that we here have considered a nonrotating atmosphere; otherwise dfi^dt 
and di /dt would be nonzero . We have also bypassed the complexities of assuming an exponential 
density profile by averaging the effects over one revolution of the satellite. 


OBLATENESS 


In this section we will merely list the variational equations due to the oblateness of the central 
body. The equaticms in and are derived in detail in Reference 11 and those in Jg and were 
taken from Reference 12. 


Equations in J^, J/ 



a 0 


45n J 2 
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(l-e^) 


siii^ i sin 2u/ 


1^ 2-\ 

[is " M 


3n JjR,' 


9nJ 2 r4 


2 *^6 

- cos i cos i 

2p2 4P^ 


i [t ‘ i (f + "Jr 0 


(91) 

(92) 


/ 7 5 . 2 \ _ 

\T5* “ T" ^ j ® cos 2u) 


( 93 ) 
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. 75 . . A 

1 7 - 55 in"' 1 f sin'» ij 


cos 2u.) 


^7 - *y‘ sin*'^ sin^ i + c*^ ^7+5 


- 2 • '♦S . 4 

sxrr' 1 “ sin^ ij 


(95) 


where 

is the equatorial radius of the planet 

p = a (i - e^) is the semilatus rectum and the subscripts o, Jj mean oblateness and Jj terms 

only. 

Equations in J , 


where 


50.J3 


2P tan i 




o.j, 


\ n R^^ J, 


cos sin i (5 cos^ i ” ^) 


o.j, 


3nR3 j 


a J3 


8P3 


e s in 6j cot i (l5 cos ^ i - 1 1) 


o.j, 


3nR,^ J3 
8P3 


e cos o> cos i (S cos ^ i - l) 


'O.J, 


3oRg J 3 +4e^) , . . a . ,\ /dfl') 

__ ;in i (5 cos^ 1 - 1 j ~ cos 1 

0.J3 


8P^ 


sino) sin 



(96) 

(97) 

(98) 

(99) 

( 100 ) 

( 101 ) 
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Equations An J4 


/ de\ 

/ di\ 

'‘‘'•/o.J, 


2P tan i 


, /dl\ 
^ \dt/ 


0‘J4 


ISnR,-' Ja ^ , 

^ 2 tx sin^ i (7 cos^ i *• 1 ) 


15 nR;> J, 

aap^ 

ISnR; J< 
64P^ 


cos 


i ^2 (7 cos^ i - 3 ) + ^7 cos^ i - I ♦ 4 8in^a>(7 cos^ i ~4)j^ 


sin 2wsin 2i (7 cos^ i " i) 


1^1 

ISnR'J, f 


16P^ \ 


(102) 


(103) 


(104) 


(105) 


|s “ 28 sin^ i + 21 sin^ i - sin^ o> sin^ i (7 cos^ i - l) 

+ ©2 jg ^ sin^ i + sin^ i + sin® ~ 35 sin^ i + ~y sin^ i^ \ 


(106) 


where a definition identical to Equation 101 holds with the proper change in the subscripts. 


THE TOTAL VARIATIONAL EQUATIONS 

As shown in Reference 3, if a ^ubly-averaged solar disturbing function is used, then the re- 
sulting set of equations separates into two uncoupled systems of 3 equations each. This set of un- 
coupled equations was used in Reference 13 along with the variation in eccentricity obtained with a 
singly- averaged disturbing function. Howevqr that study included only tiiird-body effects. Using 
these uncoupled equations, one can obtain initial conditions for various types of behavior due to 
third-body effects. Williams and Lorell in Reference 3 give a mathematical descripticm of the 
motion and Uphoff in Reference 4 gives a very good description of the interac^ns of the varicnis 
perturlmtions. Using the doubly-averaged third-body disturbing funicticm and superimposing medium 
period terms, Uphoff defines the regions in which Hie doubly- averaged equations yield very accurate 
results and thmi defines those regicms where one must resort to the singly-averaged efquaticms. A 
simple analysis of Icmg-term third-body effects may also be obtained by merely looking at these 
doubly- averaged terms. A value of 135® or 315® in w will result in long-period third-body perturba- 
tions that will raise periapsis; cj = 45® or 225® will lower periapsis and no change is obtained for 
^ = 0°, 90° or 180°. As Uphoff points out, the Icxng periochc variatiims in eccentricity are more 
prominent at high inclinaticms, while at lc}w inclinations the medium -period terms are dominant 
Also the medium-period variations increase as the eccentricity grows under the influence of the 



long-period perturbations. It is unfortunate that when one goes to the singly-averaged equations 
the ability to obtain tliis type of information analytically is lost 

For the purposes of determining the total variations, the equations may now be written in the 
following manner. 


(In 

at 


cie 

It 


do 

dt 


di 

dt 


W 


dM 

It 


r /da\ 

~ ^ *' VcTt*] 

\«Vo.j 3 

/da\ 

■Wo,;, 


(107) 

/de\ T 

. . e 4 

/de\ 

Wh.i, * wic.,, * Wo.j, 


(108) 

/dQ\ (dTi 

" Vdt^ ^ 

\ 1^ 
'0.J2 ^ 


(109) 

/di\ /dO 

1.,, ^ 


( 110 ) 

/ da>\ f dki 

" \HFj 3 + [it 

\ /dcci\ 

/O.J 2 ^ \^/OJ3 V^Wo.J^ 


( 111 ) 
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In these equations one may include or not include any of the indicated terms in order to determine 
the effects of each. However with the exception of Equation 112 (which is included only for com- 
pleteness), the above system yields an accurate description of the average motion of the orbit. As 
discussed earlier one of the princif^l advantages of this system is that the interaction of the medium 
and long-period terms is treated in a more realistic manner. This becomes particularly important 
in the case of Mars, where oblateness and third- body effects can combine to cause very complex 
behavior. 

SAMPLE CASES 

Figures 3 and 4 show a comparison between the variational method described here and an Enc^e 
n-body numerical integration program. The slight differences noticed here are due to short-period 
terms wMch are present in the n-body program and to diffm'^ices in constants. The curve in Fig- 
ure 3 produced by the n-body program tocdc more than one hour of 360/95 time, while the variational 
method took 30 seconds on the same machine and only half of this was iq)ent in the computation step. 

Figure 5 from Reisrmice 4 shows a cx)nmarison with a Cowell n-body program. Again the 
agreement betweei the two is excellm^t 
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Figures 6 and 7 show the importance of including the atmosphere of a Martian orhiter, partic- 
ularly when a long lifetime is required. 

These few sample cases illustrate the usefulness of the program outlined above. The extreme 
speed and reasonable accuracy should prove the program to be a VfUuable aid in determining the 
orbital behavior for a multitude of initial conditions. 
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Figure 3— Mars; Periopsls vs time. Encke and variational comparison. 
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Figure 4— Mars* Periopsis radius vs time, Encke and variational comparison. 
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Figure 5— Mars? Ppriopsis vs time. Cowell and voriational comparison. 
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APPENDIX A 

Mars and Venus Transformations 


Mars Transformation 

To calculate Mars-centered position of the Sun at time t in mean Mars equator of date (t). 
Given the Earth- centered position in mean equinox and Earth equator of date. 

1. Positian of Mars’ Pole at beginning of year t ^ . 

s 317.“793416667 + 0^6520833 X (t' - 1950.0) 

8 o' ^ 54J0575 + 0^85 X 10 (t' - 1950.0) 

Secular variation 

(Julian date)j - (Julian date)^» 

T = 365. 25 


then 

ao ^ ao* - 0f001013 T , 

So So' - 0^000631 T 

2. Determine H; The longitude of the ascending node of the Mars orbit along the ecliptic from 
the vernal equinox of the Earth. 

# 

(Julian date)j - 2415020.0 
“ 36525.0 

Cl ^ 48^78644167 + 0.°77099i67 T, - 0^13888889 X 10"* 

Calculate i : The inclinaUcui of th ; Martian orMt plane to the ecliptic plane. 

„ ' ' i, i 

i lf85033S338 - 0^«75 x lO'^ t„ + 0.12611111 X 10--» T* 

6 e 

4. Calculate e ; The obliquity of the ecliptic. 

£ = 2S°45229444 - 0.°130125 x 10' » T. - 0“16S88889 x 10 T ^ + 0f50277778 x 10‘® Tj» 
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5. Calculate Unit Vector o along pole of Mars 




08 'Xq cos 


0* 


sin U(j cos 


Siqi Sin 


6. Calculate Unit Vector Q perpendicular to orbit of Mars 

b - (sin i sin \l, - sin i cos H cos e - cos i sin t » - sin i cos 0 sin c + cos i cos c) 

Calculate 6 x e unit vector to ascending node of equator of Mars on the orbit of Mars 
(autumnal equinox). 

Calculate eq unit vector to the ascending node of the equator of Mars on the equator of 
Earth. 

EQ = (“sinajj, cosa^, o) 

9. Calculate ioy . The arc of the equator of Mars from its ascending node on the equator of 
Earth to its ascending node on the orbit of Mars. 

cos = EQ • (Ox e) 
sinw^ ~ |EQx(Oxe)l 

10. Calculate longitude of the ascending node of the equator of Mars on the equator of the 
Earth measured from Earth’s vernal equinox. 

' <Xq + rr/2 ^ 

11. Calculate i^ ; The inclination of the Mai*tia^ equator with respect to the Earth's equator. 

irf = ^/2 - Sg 

12. Calculate Rotation Matrix (P) 

' ' """ ■' |i 

Pjj = cosoj^cosH^ - sin6j^ sinH^ cos i^ 

Pj 2 = coscu^sinn^+ sinw^'cosn^ cos 

I! 

Pj j = sin (w^sin i^ \ 

0 

Pjj = ■" sinoJ^cosH^ - cos ( 1 ;^^ sin (2^ cos i 
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Pj 2 = “ f in f in + cog coa cot i ^ 

Pjj = cos sin i j 
P31 “ ainn^sini^ 

P32 = - cos n^sin 

P33 ^ cosi^. 

Let be a vector from the Earth to Msurs in mean Earth equator and equinox date (t ). Let 
be a vector from the Earth to the Sun in mean Earth equator and equinox of date. Itai 

f = ^«p “ %cf • 

Let be the position of the fikm in Mars mean equator and equinox of date. Then 

7 ^ (P)r 

also 

Vrf = (P)V . 


X axis is now towards ascending node of equator of Mars on orbit of Mars. 
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Venus Trans formaticm (Taken frcwn Reference 14). 


Position of Pole 

a = 98 - 0.0015551 (t - 1964.5) deg. 

8 « -69 - 0.0007748 (t - 1964.5) deg. 
t - 1964.5 in tropical years 
A s 180f075 (1964.5) 


Let A be rotation matrix from mean Earth equator and equinox 1952 to mean Earth equator and 
equinox of date. Then 


X = S3 Sj Sj AX^ 950 0 


where x is Venus equator and equinox of date 




